
########################################################################
#      File Name:	LRSS table appendix.R							                   #
#      Date:   	June 2025     			            				               #
#      Author: 	Juan Acevedo-Ossa		                        	         #
#      Purpose:   Replication file for "The Indirect Effects of	       #
#				          Structural Power										                 #
#      Software: RStudio 2024.09.1                                     #
#      Input Files:  GVC.csv, panel_unga.csv"                      	   #	
#      Output File: "LRSS_table.docx"                                  #
#      Machine:  16.0 GB RAM - Intel Core i7 @ 2.50 GHz 16 cores       #
########################################################################



rm(list = ls())
#---------------------------#

#---------------------------#
# Load required packages
#---------------------------#


#install.packages("foreign")
#install.packages("ggplot2")
#install.packages("MASS") 
#install.packages("dplyr")
#install.packages("texreg")
#install.packages("sf")         
#install.packages("spdep")      
#install.packages("spatialreg") 
#install.packages("flextable")
#install.packages("officer")



library(foreign)
library(ggplot2)
library(MASS) 
library(dplyr)
library(texreg)

library(sf)         
library(spdep)      
library(spatialreg) 
library(flextable)
library(officer)

#seed for simulation
set.seed(12345)

#---------------------------#




# ----------------------------------- #
# Load data and setup
# ----------------------------------- #

# Get the directory where the script is located
path <- "C:/Users/juanp/OneDrive - City University of New York/The Indirect Effects of Structural Power"
setwd(path)


unga_panel <- read.csv("./Data/panel_unga.csv")
wp <- as.matrix(read.csv("./Data/GVC matrices/GVC.csv"))
np <- nrow(unga_panel) # observations in panel


# Generate unique panel-id name (countrycodes_year) to use in neighbors and weights
unga_panel$id <- paste(unga_panel$countrycodes, unga_panel$year, sep = "_")
head(unga_panel$id)

# Create a neighbor listw object for pn data
# Note - "nbs_lw_pn" = neighbors_listw_panel
nbs_lw_pn <- mat2listw(x = wp,
                       row.names = unga_panel$id,
                       style     = "W")

# Update w to row-standardized [relevant for marginal effect calculations!]
wp <- listw2mat(listw = nbs_lw_pn)
colnames(wp) <- rownames(wp)

# For LRSS estimates:
n  <- 152                 # Cross-section units
# balanced panel - 152 observations per year
w  <- wp[1:n, 1:n]  

# ----------------------------------- #
##  Spatiotemporal lag (STAR) model  
# ----------------------------------- #

unga.lagsartslm <- lagsarlm(formula = alignusimp  ~ lalignusimp + sanct_stock  + ltradegdpuslog  + laidusgdplog  + polyarchy  + pop_log  + loggdppercapita + corruption_control  + imf  + atop + right_wing + dip_visits_usa + unsc_non_p + developing,
                            data    = unga_panel,
                            listw   = nbs_lw_pn)
summary(unga.lagsartslm)

screenreg(unga.lagsartslm)


# ----------------------------------- #
##       LRSS - Estimation          ----
# ----------------------------------- #

betas = c('sanct_stock', 'ltradegdpuslog', 'laidusgdplog', 'polyarchy', 'pop_log', 
          'loggdppercapita', 'corruption_control', 'imf', 'atop', 'right_wing', 'dip_visits_usa', 
          'unsc_non_p', 'developing')

Id_cs <- diag(1,n,n) 
rho   <- coef(unga.lagsartslm)['rho'] 
phi   <- coef(unga.lagsartslm)['lalignusimp'] 

LRSSEffects <- list()

for (bt in betas){

  beta  <- coef(unga.lagsartslm)[bt]  
  Id_cs <- diag(1,n,n)  
  
  #LRSS formula
  M_lrss  <- solve(Id_cs - rho*w - phi*Id_cs)
  bM_lrss <- beta * M_lrss
  
  # Spatial Effects:
  dir.eff <- mean(diag(bM_lrss))
  ind.eff <- (sum(bM_lrss) - sum(diag(bM_lrss))) / n
  tot.eff <- sum(bM_lrss) / n
  vcvm <- vcov(unga.lagsartslm)[c('rho','lalignusimp',bt),
                                c('rho','lalignusimp',bt)]
  
  #LRSS CIs - Delta Method
  
  bZWZ <- beta * w %*% M_lrss  # d_cf_rho1
  bZZ  <- beta * M_lrss        # d_cf_phi1
  
  # Create storage matrix for variance estimates:
  delta <- matrix(0, nrow = n, ncol = n)
  
  for (i in 1:n){
    d_cf_rho  <- bZWZ[,i]                           # Gradient w.r.t. rho
    d_cf_phi  <- bZZ[,i]                            # Gradient w.r.t. phi
    d_cf_beta <- M_lrss[,i]                         # Gradient w.r.t. beta
    d_cf      <- cbind(d_cf_rho,d_cf_phi,d_cf_beta) # Collect gradients in NX3 matrix
    delta1    <- d_cf %*% vcvm %*% t(d_cf)          # Delta Method vcov(LRSS) Estimate
    delta[,i] <- diag(delta1)                       # Diagonal of which are DM Variances
  }
  
  DMse     <- sqrt(delta)  # Convert   NxN variance estimates to standard errors
  DMtstats <- bM_lrss/DMse      # Calculate NxN matrix of t-statistics
  
  # Update row/column names for delta method standard error, tstat matrices
  colnames(DMse)     <- unga_panel$countrycodes[1:n]
  rownames(DMse)     <- unga_panel$countrycodes[1:n]
  colnames(DMtstats) <- unga_panel$countrycodes[1:n]
  rownames(DMtstats) <- unga_panel$countrycodes[1:n]
  
  DMse[1:5,1:5]
  DMtstats[1:5,1:5]
  
  
  #-------------------------------CIs LRSS
  
  sims <- 1000   # Set the number of iterations to simulate
  
  # Create array to store n*n matrices over 1000 sims
  sim.effs <- array(data     = 0,
                    dim      = c(n,n,sims),
                    dimnames = list(unga_panel$countrycodes[1:n],
                                    unga_panel$countrycodes[1:n]))
  
  # Simulate
  for(i in 1:sims){
    draws <- mvrnorm(n         = 1,
                     mu        = c(rho, phi, beta),
                     Sigma     = vcvm)
    
    sim.rho  <- draws['rho']
    sim.beta <- draws[bt]
    sim.phi  <- draws['lalignusimp']
    
    sim.m  <- solve(Id_cs - sim.rho * w - sim.phi * Id_cs)
    sim.bm <- sim.beta * sim.m
    
    rownames(sim.bm) <- unga_panel$countrycodes[1:n]
    colnames(sim.bm) <- unga_panel$countrycodes[1:n]
    
    sim.effs[,,i] <- sim.bm
  }
  
  ade.median <- median(apply(sim.effs, 3, function(x) median(diag(x))))
  ade.ci   <-   quantile(apply(sim.effs, 3, function(x) median(diag(x))),
                         probs = c(0.025, 0.975))
  
  LRSSEffects[[bt]] <- cbind('Variable' = bt,
        'Median Response' = ade.median,
        'Lower Bound' = as.numeric(ade.ci[1]),
        'Upper Bound' = as.numeric(ade.ci[2]))
  
}

# For table of LRSS effects
results_df <- data.frame(Variable = character(), 
                         Median_Response = character(), 
                         Lower_Bound = character(), 
                         Upper_Bound = character(), 
                         stringsAsFactors = FALSE)

for (item in LRSSEffects) {
  results_df <- rbind(results_df, data.frame(
    Variable = item[1],
    Median_Response = item[2],
    Lower_Bound = item[3],
    Upper_Bound = item[4],
    stringsAsFactors = FALSE
  ))
  
}


print(results_df)
# Create a flextable object from the data frame
ft <- flextable(results_df)

# Format the flextable (optional)
ft <- set_header_labels(ft, 
                        Variable = "Variable", 
                        Median_Response = "Median Response", 
                        Lower_Bound = "Lower Bound", 
                        Upper_Bound = "Upper Bound")

ft <- autofit(ft)

# Export the flextable to a Word document
doc <- read_docx() %>%
  body_add_flextable(value = ft) %>%
  body_add_par(" ", style = "Normal") # Add a space after the table

# Save the Word document
print(doc, target = "LRSS_table.docx")
